import math
import gamma
from matplotlib import pyplot

def loadSizeDistribution():
    name = '/home/jonk/Documentos/papers/Bombard/HQIronParticleSizeDistribution.dat'
    size = [0.0]
    freq = [0.0]
    for line in file(name):
        line = line.replace('\n','')
        val=line.split(' ')
        res=[]
        for v in val:
            if v!='':
                res.append(float(v))
        size.append(res[0])
        freq.append(res[1])
    #Normalizacion
    suma=0.0
    for i in range(1,len(freq)):
        suma += (size[i]-size[i-1])*(freq[i]+freq[i-1])*0.5
    freqnorm=[]
    for f in freq:
        freqnorm.append(f/suma)
    return [size,freqnorm]
                
    
def schultz(x,xm,alpha):
    C1=(alpha+1.0)/xm
    uno = math.pow(C1,alpha+1)/gamma.Gamma(alpha+1)
    fun = uno*math.pow(x,alpha)*math.exp(-C1*x)
    return fun

xm=0.8
alpha =0.4
x=[]
y=[]
desde=0.0
hasta=10.0
puntos=100
paso =(hasta-desde)/float(puntos)
for i in range(puntos):
    xaux=float(i)*paso
    x.append(xaux+0.3)
    y.append(schultz(xaux,xm,alpha))

pyplot.plot(x,y,'r-')
#pyplot.show()

[size,freq] = loadSizeDistribution()
pyplot.plot(size,freq,'ro')
pyplot.show()